import matplotlib.pyplot as plt
import networkx as nx
from scipy import io
import pandas as pd
import numpy as np
import scipy
import json
from cvxpy import *
from lib import *
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
A=scipy.io.loadmat("./A1.mat")["A1"]
xy=scipy.io.loadmat("./xy.mat")["xy"]
y=scipy.io.loadmat("./y.mat")["y"]
To perfectly reconstruct a signal x from a sample r
given a sample, we want to solve
where the last equation is valid because the signal is band limited. So, to perfect recovery (, we need
and this inverse exists if and only if the recovery condition is satisfied and the bigger singular value of the matrix is smaller than one.
L, U, lambdas, Px, Py, y, graph_A = preprocess(A, xy, y)
plt.figure(figsize=(60,40))
c_min=min(y)
c_max=max(y)
nx.draw_networkx(graph_A, pos=xy, arrows=True, with_labels=False, edge_color="green", node_color=y, node_size=400, node_shape="o", font_size=30,
linewidths =10, width=2, cmap=plt.cm.jet, vmin=c_min, vmax=c_max)
cmap=plt.cm.jet
vmin = c_min
vmax = c_max
sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array(y)
plt.colorbar(sm)
plt.grid()
plt.title("traffic flow signal", fontsize=100)
plt.show()
plt.figure(figsize=(16,6))
plt.plot(y, linewidth=3)
plt.title("traffic flow", fontsize=30)
plt.xlabel("nodes", fontsize=20)
plt.ylabel("flow", fontsize=20)
plt.grid()
plt.figure(figsize=(16,6))
plt.hist(y, bins=20, edgecolor='white')
plt.ylabel("frequency")
plt.xlabel("flow")
plt.grid()
plt.figure(figsize=(16,6))
plt.plot(lambdas, "o")
plt.ylim([0,0.2])
plt.axhline(lambdas[14], color="red")
plt.axhline(lambdas[15], color="red")
plt.axvline(14, color="green")
plt.xlim([0,20])
plt.grid()
lambdas[0]
U[:,0]
clusters, clusters2 = clustering_generator(U)
plt.figure(figsize=(50,40))
nx.draw_networkx(graph_A, pos=xy, arrows=True, with_labels=False, edge_color="green", node_color=clusters, node_size=400, node_shape="o", font_size=30,
linewidths =10, width=2)
plt.title("traffic flow", fontsize=30)
plt.grid()
plt.show()
plt.figure(figsize=(50,40))
nx.draw_networkx(graph_A, pos=xy, arrows=True, with_labels=False, edge_color="green", node_color=clusters2, node_size=400, node_shape="o", font_size=30,
linewidths =10, width=2)
plt.title("traffic flow", fontsize=30)
plt.grid()
plt.show()
plt.figure(figsize=(50,40))
u=U[:,14]
c_min=min(u)
c_max=max(u)
nx.draw_networkx(graph_A, pos=xy, arrows=True, with_labels=False, edge_color="green", node_color=u,
node_size=400, node_shape="o", font_size=30,
linewidths =10, width=2, cmap=plt.cm.jet, vmin=c_min, vmax=c_max)
cmap=plt.cm.jet
vmin = c_min
vmax = c_max
sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array(u)
plt.colorbar(sm)
plt.title("traffic flow", fontsize=30)
plt.grid()
plt.show()
L_approx = np.dot(np.dot(U, np.diag(lambdas)), U.T)
L_approx[np.abs(L_approx) < np.exp(-16)] = 0
(L_approx - L).sum()
bandwidth = [10, 20, 30, 40, 50]
y_f = np.dot(U.conj().T, y)
y_f.shape
plt.figure(figsize=(16,6))
plt.plot(y_f, linewidth=2)
plt.title("graph fourier transform")
plt.ylabel("gft")
plt.xlabel("frequency")
plt.grid()
y_f_sorted = [(i,np.abs(j)) for i,j in enumerate(y_f)]
y_f_sorted = sorted(y_f_sorted, key=lambda u:u[1], reverse=True)
F_global = [i[0] for i in y_f_sorted]
for b in bandwidth:
plt.figure(figsize=(16,6))
plt.plot(y_f)
y_f_sample = np.zeros(len(y_f)) * np.NaN
f_list = F_global[:b]
y_f_sample[f_list] = y_f[f_list]
plt.plot(y_f_sample, "o")
plt.title("bandwidth: " + str(b), fontsize=20)
plt.grid()
v = np.zeros(len(y)) #np.random.normal(0, 0.1, len(y))
v = v.reshape(len(v), 1)
choose the best sampling set is a discrete combinatorial problem: an exaustive search is infeasible. To deal with the exponential complexity of this problem we have two possibilities:
in this work we concentrate our attention on different greedy strategy for local optimal solution.
def random_routine(U, bandwidth, y, F_global):
Random = {b:[] for b in bandwidth}
a = [i for i in range(len(y))]
for b in bandwidth:
F = b
S = F + 3
F_list = F_global[:F]
Uf = U[:, F_list]
index_list = np.random.choice(a, size=S, replace=False)
print(index_list)
y_reconstructed, y_sample, matrix_rec, y_sampled = y_reconstruction_routine(U, y, index_list, F_list)
y_reconstructed[y_reconstructed < 0] = 0
error = compute_error(y_reconstructed, y)
Random[b] = [index_list, y_reconstructed, error, y_sample, matrix_rec, y_sampled]
print(b)
return Random
Random=random_routine(U, bandwidth, y, F_global)
index=["index_list", "y_reconstructed", "error", "y_sample", "matrix_rec", "y_sampled"]
Random_df=pd.DataFrame(Random, index=index)
Random_df.loc["error"]
def routine_without_noise(U, bandwidth, y, F_global):
res = {b:[] for b in bandwidth}
for b in bandwidth:
F = b
S = F + 3
F_list = F_global[:F]
Uf = U[:, F_list]
index_list = greedy_e_opt(Uf, y, S)
try:
y_reconstructed, y_sample, matrix_rec, y_sampled = y_reconstruction_routine(U, y, index_list, F_list)
print(index_list)
y_reconstructed[y_reconstructed < 0] = 0
error = compute_error(y_reconstructed, y)
except np.linalg.LinAlgError:
print("error")
error=np.NaN
y_reconstructed_e=np.NaN
y_sample=np.NaN
matrix_rec=np.NaN
y_sampled=np.NaN
res[b] = [index_list, y_reconstructed, error, y_sample, matrix_rec, y_sampled]
print(b)
return res
E_opt=routine_without_noise(U, bandwidth, y, F_global)
E_opt_df = pd.DataFrame(E_opt, index=index)
E_opt_df.loc["error"]
def routine_with_noise(U, bandwidth, y, F_global, v, t="e"):
res_noise = {b:[] for b in bandwidth}
for b in bandwidth:
F = b
S = F + 3
F_list = F_global[:F]
Uf = U[:, F_list]
if t == "e":
index_list = greedy_e_opt(Uf, y, S)
elif t == "a":
index_list = greedy_a_opt(Uf, y, S, v)
elif t == "d":
index_list = greedy_d_opt(Uf, y, S, v)
try:
y_reconstructed, y_sample, matrix_rec, y_sampled = y_reconstruction_routine_noise(U, y, index_list, F_list, v)
y_reconstructed[y_reconstructed < 0] = 0
print(index_list)
error = compute_error(y_reconstructed, y)
except np.linalg.LinAlgError:
print("error")
error=np.NaN
y_reconstructed=np.NaN
y_sample=np.NaN
matrix_rec=np.NaN
y_sampled=np.NaN
res_noise[b] = [index_list, y_reconstructed, error, y_sample, matrix_rec, y_sampled]
print(b)
return res_noise
E_opt_noise=routine_with_noise(U, bandwidth, y, F_global, v, "e")
A_opt_noise=routine_with_noise(U, bandwidth, y, F_global, v, "a")
D_opt_noise=routine_with_noise(U, bandwidth, y, F_global, v, "d")
E_opt_noise_df=pd.DataFrame(E_opt_noise, index=index)
A_opt_noise_df=pd.DataFrame(A_opt_noise, index=index)
D_opt_noise_df=pd.DataFrame(D_opt_noise, index=index)
E_opt_noise_df.loc["error"]
A_opt_noise_df.loc["error"]
D_opt_noise_df.loc["error"]
plt.figure(figsize=(20,10))
#r = [Random[b][2] for b in bandwidth]
e = [E_opt_noise[b][2] for b in bandwidth]
a = [A_opt_noise[b][2] for b in bandwidth]
d = [D_opt_noise[b][2] for b in bandwidth]
#plt.plot(bandwidth, r)
plt.plot(bandwidth, e, linewidth=4)
plt.plot(bandwidth, a, linewidth=4)
plt.plot(bandwidth, d, linewidth=4)
plt.legend(["E_opt", "A_opt", "D_opt"], fontsize=20)
plt.title("error vs bandwidth", fontsize=20)
plt.xlabel("bandwidth", fontsize=20)
plt.ylabel("error", fontsize=20)
#plt.ylim([-10,0])
plt.grid()
#plt.savefig("error_bandwidth.pdf", format="pdf")
given a noisy sample
has to be not null for every given : this means full rank columns for and the necessary condition
and using the reconstruction formula we can obtain the reconstructed signal
def sampling_reconstruction(opt_noise, y, bandwidth, t="E"):
for b in bandwidth:
plt.figure(figsize=(16,6))
y_reconstructed=opt_noise[b][1]
y_sample=opt_noise[b][3]
error=opt_noise[b][2]
plt.grid()
plt.plot(y, linewidth=2, color="green")
plt.plot(y_reconstructed, color="red", linewidth=4)
plt.plot(y_sample, "o", color="blue")
plt.ylim([-1, 15])
plt.legend(["original", "reconstructed", "sample"], fontsize=10)
plt.title(str(t) + "_opt_noise | bandwidth " + str(b) + " | samples " + str(b+3) + " | nmse " + str(error), fontsize=20)
sampling_reconstruction(E_opt_noise, y, bandwidth, "E")
sampling_reconstruction(A_opt_noise, y, bandwidth, "A")
sampling_reconstruction(D_opt_noise, y, bandwidth, "D")
def sampling_visualization(E_opt_noise, A_opt_noise, D_opt_noise, b, title):
plt.figure(figsize=(16,6))
yE_reconstructed=E_opt_noise[b][1]
yE_sample = E_opt_noise[b][3]
yD_reconstructed=D_opt_noise[b][1]
yD_sample = D_opt_noise[b][3]
yA_reconstructed=A_opt_noise[b][1]
yA_sample = A_opt_noise[b][3]
plt.plot(yA_reconstructed)
plt.plot(yE_reconstructed)
plt.plot(yD_reconstructed)
plt.plot(yA_sample, "o")
plt.plot(yE_sample, "o")
plt.plot(yD_sample, "o")
plt.legend(["A_rec","E_rec","D_rec","a", "e", "d"])
plt.ylim([-1, 15])
plt.title(title, fontsize=20)
plt.grid()
for b in bandwidth:
sampling_visualization(E_opt_noise, A_opt_noise, D_opt_noise, b, "bandwidth " + str(b))
def sampling_graph_visualization(opt_noise, y, bandwidth, t="E"):
for b in bandwidth:
plt.figure(figsize=(20,10))
index_list = opt_noise[b][0]
y_sample = y[index_list]
xx=Px[index_list]
yy=Py[index_list]
plt.scatter(Px, Py, s=y*100, marker="o", color="green")
plt.scatter(xx, yy, s=y_sample*100, marker="o", color="red")
plt.grid()
plt.title(str(t) + "_opt_noise | bandwidth " + str(b) + " | samples " + str(b+3), fontsize=20)
plt.show()
#sampling_graph_visualization(E_opt_noise, y, bandwidth, t="E")
fig, axs = plt.subplots(2*len(bandwidth), 3, figsize=(40,100))
fig.subplots_adjust(hspace = .2, wspace=.05)
for b in range(0, 2*len(bandwidth), 2):
band = bandwidth[b//2]
y_reconstructed=E_opt_noise[band][1]
y_sample=E_opt_noise[band][3]
error=E_opt_noise[band][2]
axs[b, 0].plot(y, linewidth=2, color="green")
axs[b, 0].plot(y_reconstructed, color="red")
axs[b, 0].plot(y_sample, "o", color="blue")
axs[b, 0].set_title("E optimality| bandwidth " + str(band), fontsize=20)
axs[b, 0].set_ylim([-1, 15])
axs[b, 0].grid()
index_list = E_opt_noise[band][0]
y_sample = y[index_list]
xx=Px[index_list]
yy=Py[index_list]
axs[b+1, 0].scatter(Px, Py, s=y*100, marker="o", color="green")
axs[b+1, 0].scatter(xx, yy, s=y_sample*100, marker="o", color="red")
axs[b+1, 0].set_title("E optimality| bandwidth " + str(band), fontsize=20)
axs[b+1, 0].grid()
y_reconstructed=A_opt_noise[band][1]
y_sample=A_opt_noise[band][3]
error=A_opt_noise[band][2]
axs[b, 1].plot(y, linewidth=2, color="green")
axs[b, 1].plot(y_reconstructed, color="red")
axs[b, 1].plot(y_sample, "o", color="blue")
axs[b, 1].set_title("A optimality| bandwidth " + str(band), fontsize=20)
axs[b, 1].set_ylim([-1, 15])
axs[b, 1].grid()
index_list = A_opt_noise[band][0]
y_sample = y[index_list]
xx=Px[index_list]
yy=Py[index_list]
axs[b+1, 1].scatter(Px, Py, s=y*100, marker="o", color="green")
axs[b+1, 1].scatter(xx, yy, s=y_sample*100, marker="o", color="red")
axs[b+1, 1].set_title("A optimality| bandwidth " + str(band), fontsize=20)
axs[b+1, 1].grid()
y_reconstructed=D_opt_noise[band][1]
y_sample=D_opt_noise[band][3]
error=D_opt_noise[band][2]
axs[b, 2].plot(y, linewidth=2, color="green")
axs[b, 2].plot(y_reconstructed, color="red")
axs[b, 2].plot(y_sample, "o", color="blue")
axs[b, 2].set_title("D optimality| bandwidth " + str(band), fontsize=20)
axs[b, 2].set_ylim([-1, 15])
axs[b, 2].grid()
index_list = D_opt_noise[band][0]
y_sample = y[index_list]
xx=Px[index_list]
yy=Py[index_list]
axs[b+1, 2].scatter(Px, Py, s=y*100, marker="o", color="green")
axs[b+1, 2].scatter(xx, yy, s=y_sample*100, marker="o", color="red")
axs[b+1, 2].set_title("D optimality| bandwidth " + str(band), fontsize=20)
axs[b+1, 2].grid()
def solver(A, y):
n, m= A.shape
x = Variable(m)
b = y[:]
objective = Minimize(sum_squares(A*x - b))
constraints = [x >= 0]
prob = Problem(objective, constraints)
result = prob.solve()
x = np.array(x.value)
new = np.dot(A, x)
new = [i[0] for i in new]
return new
for b in bandwidth:
e = E_opt_noise[b][2]
Ae = E_opt_noise[b][4]
new = solver(Ae, y)
new_e = compute_error(new, y)
print("E_optimality","bandwidth: " + str(b), "greedy_error: " + str(e), "relaxed_error: " + str(new_e))
a = A_opt_noise[b][2]
Aa = A_opt_noise[b][4]
new = solver(Aa, y)
new_a = compute_error(new, y)
print("A_optimality","bandwidth: " + str(b), "greedy_error: " + str(a), "relaxed_error: " + str(new_a))
d = D_opt_noise[b][2]
Ad = D_opt_noise[b][4]
new = solver(Ad, y)
new_d = compute_error(new, y)
print("D_optimality","bandwidth: " + str(b), "greedy_error: " + str(d), "relaxed_error: " + str(new_d))
print()
we solve a convex constrained minimization problem
objective function :
def objective_function(x, Uf, Uft, Rv_inv, flag):
if flag == "d":
matrix = -log_det(Uft * diag(x) * Rv_inv * Uf)
if flag == "a":
matrix = trace(inv(Uft * diag(x) * Rv_inv * Uf)) ### inv?
if flag == "e":
matrix = -sqrt(sigma_min(diag(x) * Uf)) ### sigma_min?
return matrix
def optimization_routine(U, m, B, F_list, v, flag):
Uf = U[:, F_list]
Uft = Uf.T
Rv_inv = np.diag(np.ones(len(v)))
S = B + 3
x = Variable(m)
obj = objective_function(x, Uf, Uft, Rv_inv, flag)
objective = Minimize(obj)
constraints = [sum_entries(x) == S, 0 <= x, x <= 1]
prob = Problem(objective, constraints)
result = prob.solve()
x = np.array(x.value)
return x
def index_selection(x, B):
S=B+3
new = [ (i[0], j) for j, i in enumerate(x)]
new = sorted(new, reverse=True)
index_list = [i[1] for i in new[:S]]
return index_list
def relaxed_solver(U, y, m, B, F_list, v, flag):
x = optimization_routine(U, m, B, F_list, v, flag)
index_list = index_selection(x, B)
y_r, y_sample, matrix_rec, y_sampled = y_reconstruction_routine_noise(U, y, index_list, F_list, v)
return y_r
for b in bandwidth:
m = len(y)
F_list = F_global[:b]
# e = E_opt_noise[b][2]
# new = relaxed_solver(U, y, m, B, F_list, v, "d")
# new_e = compute_error(new, y)
# print("E_optimality","bandwidth: " + str(b), "greedy_error: " + str(e), "relaxed_error: " + str(new_e))
# a = A_opt_noise[b][2]
# new = relaxed_solver(U, y, m, B, F_list, v, "d")
# new_a = compute_error(new, y)
# print("A_optimality","bandwidth: " + str(b), "greedy_error: " + str(a), "relaxed_error: " + str(new_a))
d = D_opt_noise[b][2]
new = relaxed_solver(U, y, m, B, F_list, v, "d")
new_d = compute_error(new, y)
print("D_optimality","bandwidth: " + str(b), "greedy_error: " + str(d), "relaxed_error: " + str(new_d))
print()